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(N 1 ABSTRACT 

The non-linear fluid dynamics of a warped accretion disc was investigated in an earlier 
paper by developing a theory of fully non-linear bending waves in a thin, viscous 
disc. That analysis is here extended to take proper account of thermal and radiative 
effects by solving an energy equation that includes viscous dissipation and radiative 
transport. The problem is reduced to simple one-dimensional evolutionary equations 
for mass and angular momentum, expressed in physical units and suitable for direct 
application. This result constitutes a logical generalization of the alpha theory of 
J> ' Shakura & Sunyaev to the case of a time-dependent warped accretion disc. The local 

thermal-viscous stability of such a disc is also investigated. 
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O . 1 INTRODUCTION 

■ There are many good reasons for investigating the dynamics of thin accretion discs that are warped or twisted, in the sense 
Or that fluid elements closely follow circular Keplerian orbits, but the orbital plane varies slowly with radius and possibly with 
Q ■ time. The existence of such discs is supported by a growing body of observational evidence and theoretical reasoning. 

Direct observational evidence for warped profiles has been obtained in exceptional casesQ, including the circumnuclear 
C/j ■ disc in the galaxy NGC 4258, revealed through very-long-baseline interferometric observations of water masers (Miyoshi et al. 
. . \ 1995), and the dusty circumstellar disc of (5 Pic (e.g. Heap et al. 2000). Indirect observational evidence comes from several 
^ ■ low-mass X-ray binaries in which a long periodicity or other modulation is attributed to a precessing tilted disc (see Wijers & 
Pringle 1999 and references therein). The best example is the eclipsing system Her X-l (Tananbaum et al. 1972; Katz 1973; 
» I ■ Roberts 1974), but the precessing jets of SS 433 (e.g. Margon 1984) also presumably originate in a tilted disc. Further indirect 
evidence comes from pre-main-sequence binaries or unresolved young stellar objects from which two jets emanate in different 
directions (see Terquem et al. 1999 and references therein). If the jets are emitted perpendicular to the discs that produce 
them, at least one of the discs must be misaligned with the binary orbital plane. Such a misalignment appears to have been 
captured directly in images of HK Tau (Stapelfeldt et al. 1998; Koresko 1998). In each of these cases, the disc cannot simply 
be rigidly tilted, but must become warped under the influence of the tidal potential. 

Even without observational motivation, there are theoretical grounds for considering warped discs. In any situation in 
which, for historical reasons, the disc is misaligned with the spin axis of a central black hole (Bardeen & Petterson 1975), the 
axis of a central magnetosphere, the orbit of a companion object (Papaloizou & Terquem 1995), or the angular momentum 
vector of matter being supplied at a large distance, the disc experiences a external torque that attempts to change its 
orientation. The response of the disc is basically gyroscopic, but, unlike a rigid body, it communicates internally through 
bending waves and 'viscous' torques. In these situations the disc is likely to become warped, at least temporarily. Even in a 
system that is initially coplanar, a warp may develop spontaneously through instabilities, which may be due to resonant tidal 
interactions (Lubow 1992; Lubow & Ogilvie 2000), wind torques (Schandl & Meyer 1994), radiation torques (Pringle 1996) 
or magnetic torques (Lai 1999). 

Historically, attempts to describe the dynamics of warped accretion discs have involved deriving linear equations for 
slowly varying m = 1 bending disturbances of an initially flat disc. It is usually assumed that the agent responsible for 

* The more familiar examples of warped spiral galaxies are not considered here, because they are likely to be controlled by different 
physical effects. 
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angular momentum transport in accretion discs (probably magnetohydrodynamic turbulence) may be treated as an isotropic 
effective viscosity, described by the dimensionless parameter a. For Keplerian discs an important qualitative distinction exists 
between significantly viscous discs (a J; H/r, where H/r is the angular semi-thickness of the disc), in which warps evolve 
diffusively (Papaloizou & Pringle 1983), and almost inviscid discs (a ^ H/r), in which warps propagate as approximately 
non- dispersive waves (Papaloizou & Lin 1995). 

However, a non-linear theory is required in many cases where the disc is warped sufficiently to have observational 
consequences. Different approaches have been taken to the non-linear problem. Direct numerical simulations have been carried 
out using smoothed particle hydrodynamics (SPH; e.g. Larwood et al. 1996; Nelson & Papaloizou 1999). This method makes 
no attempt to resolve the magnetohydrodynamic turbulence but involves an artificial viscosity over which some control can 
be exerted. It is most suitable for relatively thick discs of limited radial extent, and these investigations have focused mainly 
on the case a <y H/r. Usually good agreement is found with linear theory, but new, non-linear effects have also been found. 
The computational expense makes it difficult to follow the evolution of the disc on the long, viscous time-scale. 

Meanwhile, the semi-analytical method has been developed into a non-linear theory. Pringle (1992) adopted plausible 
evolutionary equations for a thin, viscous disc with a warp of arbitrary amplitude. The equations are simple conservation laws 
for mass and (vector) angular momentum in one (radial) dimension. The advantages of this approach are that the equations 
have a clear physical interpretation in terms of torques acting on rings, and that the numerical implementation is much less 
demanding than for a direct numerical simulation, so that the evolution of the disc on the viscous time-scale is easily described. 

In an effort to provide a more complete justification of the latter approach, I have developed a theory of fully non-linear 
bending waves in a thin, viscous disc, starting from the Navier-Stokes equation in three dimensions (Ogilvie 1999, hereafter 
Paper I). It was shown that, in cases where the resonant wave propagation found in inviscid Keplerian discs (Papaloizou & 
Lin 1995) may be neglected, the problem can indeed be reduced to one-dimensional evolutionary equations as suggested by 
Pringle (1992). Moreover, the dimensionless coefficients in the equations, and their variation with the basic parameters of the 
disc and with the amplitude of the warp, were derived, and a proper connection was made with the existing linear theory. 
Further discussion of the meaning and validity of this approach is given in Section 5.2 below. 

It was recognized that the treatment of thermal physics in Paper I was oversimplified. In order to make the non-linear 
problem tractable, a polytropic relation was assumed to hold locally in radius. In reality, the thickness of the disc and the 
stresses between neighbouring rings depend on a balance between the local rate of dissipation of energy and the radiation 
from the surfaces of the disc. In a warped disc, as a result of induced motions, the viscous dissipation is enhanced where the 
warp is strongest, and this affects the non-linear dynamics of the warp. The aim of the present paper is to give a proper 
treatment of thermal physics by solving an energy equation including viscous dissipation and radiative transport. Remarkably, 
the problem can still be reduced by separation of variables to a set of dimensionless, non-linear ordinary differential equations 
(ODEs) which yield the required coefficients. 

The remainder of this paper is organized as follows. I present the extended set of equations in Section 2 and solve for the 
instantaneous azimuthal and vertical structure of the warped disc and the induced motions. I then deduce the complete set 
of evolutionary equations for the warped disc. In Section 3 the evolutionary equations are tested for stability with respect to 
short-wavelength perturbations. The outcome of numerical evaluation of the dimensionless torque coefficients and the local 
stability of the disc is presented in Section 4. The results are collected in Section 5 in a concise form suitable for direct 
application, and the meaning and validity of the theory are discussed. In the Appendix, the possible effects of external 
irradiation are assessed. 

The detailed application of this theory to warped discs in X-ray binaries will be given in a forthcoming publication 
(Ogilvie & Dubus 2000). 



2 THE EXTENDED EQUATIONS AND THEIR SOLUTION 
2.1 The energy equation and radiative transport 

The basic equations and notation are taken to be the same as in Paper I (equations 1. 37-1. 39), except that the adiabatic 
condition (equation 1.38) is replaced by an energy equation that includes both viscous dissipation and radiative transport of 
energy within the Rosseland approximation. The governing equations are therefore the equation of mass conservation, 

%"<"■•- "» 

the equation of motion, 

p^j- = -Vp - pV$ + V- [(j,Vu + fi{Vuf] + V [(/it - §A*)V-u] , (2) 
and the energy equation, 
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(jvTi) Th = ~ (r=l) PV ' U + (VU) : ^ + m(Vu)T ] + " l^)( V ' u ) 2 ~ V ' F - ( 3 ) 
Here 

F = _I^!VT (4) 

is the radiative energy flux, where a is the Stefan-Boltzmann constant, T the temperature and kr the Rosseland mean opacity. 
The disc is assumed to be optically thick. The adiabatic exponent F(r) is assumed to be a prescribed function of radius, but 
a polytropic relation is no longer adopted. The equation of state is that of an ideal gas, 

kpT 

P= , (5) 

where k is the Boltzmann constant, /i m the mean molecular mass (assumed constant) and ttih the mass of the hydrogen atom. 
As in Paper I, the dynamic shear and bulk viscosities are assumed to be locally proportional to the pressure, so that 

(j, = ap/Q, (6) 
Mb = a b p/Q, (7) 

where the dimensionless coefficients a(r) and a b {r) are prescribed functions of radius. Finally, the opacity law is assumed to 
be of the general power-law form 

nn = C K p x T\ (8) 

where C K is a constant. This includes the important cases of Thomson scattering opacity (x = 0, y = 0) and Kramers opacity 
(x = 1, y = —7/2). The appropriate values of C K are discussed in Section 5.1 below. 



2.2 Set A equations 

In Paper I the equations of fluid dynamics were derived in warped spherical polar coordinates, which follow the principal 
warping motion of the disc, and were then reduced by means of asymptotic expansions for a thin disc. In physical terms, 
this separates the 'fast' orbital motion from the 'intermediate' oscillatory velocities induced by the warp and from the 'slow' 
velocities associated with accretion. The equations to be considered form two sets. 

Set A, which determines the 'intermediate' velocities, is a set of coupled non-linear partial differential equations (PDEs) 
in two dimensions, azimuthal and vertical. The physical content of Set A is as follows. In a fiat disc, the pressure and viscous 
stress are stratified in the vertical direction. When the disc is warped, this produces horizontal forces of odd symmetry about 
the mid-plane. These forces drive shearing horizontal epicyclic motions in the disc. The horizontal motions couple non-linearly 
with the warp itself to produce compressive vertical motions, causing the thickness of the disc to vary in azimuth. All the 
induced motions are subject to viscous damping and contribute to the energy dissipation rate. The problem defined by Set A 
is effectively instantaneous and local to a single annulus; derivatives with respect to time and radius do not appear because 
of the separation of scales in a thin disc. The solution of Set A must be obtained in full. 

Set B, which determines the 'slow' velocities, is a set of coupled linear PDEs with coefficients depending on the solution 
of Set A. Set B does not have to be solved, however, because the necessary information can be extracted through certain 
integrations. This procedure amounts to obtaining the solvability conditions for the set of linear equations, and these conditions 
are precisely the required equations determining the evolution of the warp. Physically, they involve a calculation of the 
instantaneous internal torque associated with the 'intermediate' velocities and the deformed structure of the disc, known from 
the solution of Set A. 

By replacing the adiabatic condition with the energy equation (^), equations (1.69) in Set A and (1.74) in Set B are 
modified. However, the manipulations of Set B in Paper I did not involve equation (1.74) and are therefore unaffected. It is 
sufficient to consider the modifications to Set A. 

Recall that (r, 9, <f>) denote warped coordinates, which are based on spherical polar coordinates, but the equatorial plane 
6 = 7r/2 is deformed to correspond to the warped mid-plane of the disc. For a thin disc, the scaled dimensionless vertical 
coordinate £ is used, where 9 = n/2 — e£ and e <C 1. Also (v r , vg, v^) denote the relative velocity components, which are 
proportional to the rates of change of the three coordinates following a fluid element. The shape of the warp is described 
through the unit tilt vector £(r,t), which has Euler angles (3(r,t) and 7(r, t). The orbital angular velocity is O(r) and the 
associated epicyclic frequency is n(r). A prime denotes differentiation with respect to r. 

In this paper a simplified notation is adopted in which {po,po, po, Pw, v r i, v$i, v^i} are written as {p, p, fi, /xb, v r , vg, v^}. 
No confusion can result because Set B will not be considered explicitly. Similarly, the new quantities {T, F, «r} represent 
their leading-order terms. Here F = —Fe represents the vertical radiative flux which is the only component present at leading 
order. 

The extended equations of Set A may be written as follows. The equation of mass conservation: 
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d(f> r <9C 



P = 



p dvg 
r~£K' 



The r-component of the equation of motion 
d v$ d 



pin 



d(j> r 



+ 



v r — 2pQ.(v r f > + rv r "(' cos f3) = — (/?' cos 4> + 7' sin /3 sin ( 



d ( dv 



p + 



(/* + |/*) 



1 \ 1 <9t>e 



+ (/?' cos + 7' sin [3 sin 0) 2 ( p-^r ) + sin0 — 7' sin /3 cos 0)^r- 



9C V a c 



The ^-component: 



p ^Q-^ — ~~~~g^J [ v e + rv r {(3' cos + 7' sin /3 sin 0)] — pQ. \rQ,C, + rv r (J3' sin 4> — 7' sin /3 cos< 



r <9C 



1 1 1 U9»9 



^2 + (/3'cos^ + 7'sin/3sin0) 2 — -j p — [vg + rv r (P' cos + 7' sin /3 sin 0)] 



-rf2(/3' cos + 7' sin /3 sin0)(/3' sin 4> — 7' sin /3 cos0)^. 



The 0-component: 



(«0 + r?v7' cos (3) + ^-t) r = + (/?' cos + 7' sin /3 sin 0)^ 



,_9 vg d \ 

- r <k ) 

+rQ'(/3' cos + 7' sin f3 sin 0) ^ . 



d_ 



V'Q£( V 'P + rV rl' COS/3) 



The energy equation: 

d vg d 



n- 



r d( 



p <9t>e 



cos + 7' sin /3 sin < 



+2^i I i -^ [«e + rtv (/3' cos <p + 7' sin /3 sin 0)] | 

f d 

+p < (/?' cos + 7' sin /3 sin 0) — [vg + rv T {(i' cos 4> + 1 sin /3 sin < 



+P 



rO.' + (J3 1 cos + 7' sin /3 sin (f>) — (v$ + rv r "f' cos f3) 



r d( 



(P' sin 4> — y sin /3 cos 0)rf2 



(w^ + rv r "f' cos /3) 



1 <9f r 



+ 



(/*-§/*) 



1 dvg 



1 <9F 



r d( J r d( 



The radiative flux: 
16ctT 3 1 dT 



(9) 



(10) 



(11) 



(12) 



(13) 



(14) 



3ftRp r d( 

Since only derivatives of the dependent variables with respect to (f> and ( appear in this set, their parametric dependence 
on r and t is neglected in the notation adopted in the following analysis. 
Recall that the dimensionless complex amplitude of the warp is 



ip — \tp\e x — r((3' + sin (3). 

The two combinations appearing in Set A are 

r(/3'cos0 + 7'sin/3sin0) = \tp\ cos(0 — x), 
r(P' sin0 — 7' sin /3cos (j>) — \tp\ sin(0 — x)- 

It will be more useful to introduce, instead of the dimensionless epicyclic frequency k, the quantity 

_ dlnfl 
9 dlnr 

These are related by k 2 = 4 — 2q. 



(15) 

(16) 
(17) 

(18) 
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2.3 Vertical structure of the disc 

It will be instructive to consider first the case of a flat disc. Then the above equations reduce to the ODEs 

-pr 2 n\, (19) 



dp 2, v . 



dC 

— = q 2 aprn, (20) 
dT _ 3K R prF 



dC ~ 16o-T 3 ' 

subject to the given equation of state, the opacity law, and the boundary conditions 

F = at C = (22) 
and 

p = T = at C = <s, (23) 

where the surface coordinate ( B is to be determined. These 'zero boundary conditions' assume that the disc is highly optically 
thick and neglect the atmosphere of the disc. The effects of finite optical thickness and external irradiation are estimated in 
the Appendix. 

Recast this problem in a dimensionless form by means of the transformations 

r( = (*U H , (24) 

p(C) = p.(C.)C„. (25) 

p(C) = P*«*)U P , (26) 

T(C) = T,(C)U T , (27) 

F(C) = F,(C.)U F , (28) 
where the starred variables are dimensionless and the U quantities are natural physical units defined by 

Uh = ( ^ a) l/(6+-2v) E ( a +-)/(6+«- a v) n -(5- 2w )/(6+«- 2 „) ( ^ H ) - (4 - V)/(6+X - 2v) ^ Ifc^ -!/(«+-*») ^ (2g) 

together with 

C/ p = El^ 1 , (30) 

C/ p = tfuhUp, (31) 

[/ T = Q^A^)^ (32) 

(7f = q 2 aQU H U p , (33) 
where 

prdC (34) 



— P*C*, (35) 
p., (36) 

-pi + "rr 3+y F», (37) 



-c. 

is the surface density. The dimensionless equations are 
dp* 
dCT 

dC* 
dT« 

p* = p*T«, (38) 
subject to the boundary conditions 

F,(0) = p»(Cs*) = T„(Cs.) = 0. (39) 
The definition of S also requires that 

p,(C.)dC. = 1. (40) 

Cb. 

For physical values of x and y it may be assumed that a unique solution exists, and can be obtained numerically. Let 
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I* = J P*(C*)C 2 d<, (41) 

be the corresponding dimensionless second vertical moment of the density. Then the dimensional quantity is 

1 = j pr\ 2 rdC =T m EU%. (42) 

The numerical solution for Thomson opacity yields ( Bt m 2.383 and 1, m 0.6777, while that for Kramers opacity gives 
Cs* ~ 2.543 and J* w 0.7094. 

This solution can now be generalized to describe the azimuthal and vertical structure of a warped disc. 



2.4 Solution of Set A equations 

In a warped disc, although many additional effects appear, the equations of Set A can be satisfied by functions having basically 
the same vertical structure as in the flat disc. However, the coefficients are different and the thickness of the disc also depends 
on <j>. As an intermediate stage to solving Set A, it is proposed that, instead of equations (^)-(^lj), the relations 

dp 
8F 

dc 

&T 3K R prF 



= -M<l>-X)pr 2 n\, (43) 
= fi(<t>-x)q 2 aprSl, (44) 

d( ~ 16aT3 ' (45) 

hold, subject to the same boundary conditions as above, where /i and /a are dimensionless functions to be determined. 
Physically, fi reflects the fact that the dissipation rate is enhanced, in a non-axisymmetric way, above its Keplerian value 
because of the additional shear associated with the induced motions. Similarly, /a represents changes to the effective vertical 
acceleration in the disc, caused by inertial forces associated with the induced motions. 
These relations can be solved as above to give 

K = C* lfi^-x)] 1/i6+x - 2y) [f2^-xT i3 - v)/i6+x - 2v) U H , (46) 
piO = P*(C*) [h^~x)r 1/i6+x - 2y) [f2(<t>~x)f- y)/(6+x - 2y) U p , (47) 
P(C) = P^)[h^^x)] in6+x - 2y Hh^~X)f +x - y)/{&+x - 2y) U p , (48) 

no = r.(c.) [h^-x)] 2n&+x - 2v) W-x)T ,(&+x - 2v) u T , (49) 

F(0 = F.(C.) [/ 1 (0- x )]( 8+ «- a «'V( 8+ «- a «')[/ a (0- x )]«/( 8+ -- a «')Ui 1 ., (50) 

where the starred functions are exactly the same as before. Note that the surface density E must be axisymmetric in order 
that mass be conserved around the orbit; this can be seen by vertically integrating equation (|^). Therefore the U quantities 
are independent of cf>. The disc thickness and the second moment are non-axisymmetric, however, with 

1 = 1* [h{4> - x )] 2 /( 6 +-" 2 ^) [f a ty _ x )]-a(s-i»)/(6+— a») E£/ 2 (51) 

Here, as in Paper I, the tilde denotes a non-axisymmetric quantity prior to azimuthal averaging. The azimuthal average is 
1=(1). 

Further propose that the velocities are all proportional to £ and are of the form 

«r(0,C) = fc(4>-X)m{, (52) 

« 9 (0,O = U(4>-x)rSK, (53) 
vM,Q + rvr{<l>,C)j'a*p = /*(0-x)rfiC, (54) 

as in Paper I. 

When this is substituted into equation (^|), one obtains 

(3 - - h(4>)fi{<t>) = {& + x- 2y)h{ct>)hWh{<t>)- (55) 

Similarly, from equation ( |l3| ) one obtains 

&(4>) = (T + l)H<t>)M<t>) + (T - ±)H<t>) (2" [h{<t>)\M cos<^] 2 + 2a [/ 4 (0) + / 3 (0)|V| cos^] 2 

+a {[f 4 (cf>) + f 3 (4>)\i>\ cos 0] |V| cos0 - |V| sin - / 3 (0)} 2 + a [-q + / 8 (^)|^| cos^.] 2 + a [/ 5 (0)] 2 

+ (a b -|a)[/4^)] 2 -g 2 a/i(0)). (56) 
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Finally, from equations (p"(i|)-([l2|) one has 

fs(4>) = H4>)h{<t>) + + [l + (Ob + §a)M4>)] M<PM\ coscb- a/ 2 (0)/ 3 (0)(l + \i>\ 2 cos 2 <t>) - af 2 (<t>M sin0, (57) 

f*(4>) = -fsttM cos + 2/a (<f>M sin + /4WO [/ 4 (0) +/ 3 (0)|V|cos0] + 1- [l + (Ob+ |«)/4(0)] / 2 (0) 

-ot-HW) [h{4>) + f3(<j>)\ip\ cos0] (1 + |V>| 2 cos 2 0) + a/ 2 (^)|i/>| 2 cos sin 0, (58) 
fs(4) = /4(0)/5(0)-(2-g)/3(0)-a/2(0)/5(0)(l + |^| 2 cos 2 0)+ga/ 2 (0)|^|cos0, (59) 

identical to equations (I.107)-(1.109). 
As in Paper I, define 

M4>-x)=i(<t>)/z, (60) 

so that 

f 6 {4>) = -2U{4>)M4>) (61) 

and 

(/«) = !• (62) 
Then the evolutionary equations for the warped disc have precisely the same form as in Paper I, viz. 

§ + I^(™E)=0 (63) 
at r or 

and 

^-(Er 2 fi£) + -^(rv^Q.t) = ~ (Q 1 lr 2 n 2 £ + Q 2 lr 3 n 2 ^± + Q a lr 3 n 2 £ x ^ , (64) 
at r or r or \ or or ) 

where v(r, t) is the mean radial velocity. The dimensionless torque coefficients Q 1 and Q4 = Q2 + iQ3 are again obtained from 
Qi = (/e l-qafa - fa ft + a/2/5 [VI cos 0] ) (65) 
and 

<3 4 = T^j ( e '*/6 [/a ~ i/s (/4 + /alV'l cos 0) + ia/ 2 (/4 + /sIV 1 cos 0)1^1 cos0 - ia/2/3 - ia/2|V>| sin0] ), (66) 

in which f n stands for / n (0). 

The relation between X and E is (from equation |jl| ) 

J=Q S X*£[4, (67) 

where 

Q g _ ^2/(6+o:-2y) 2(3— y)/ (6 + x-2y) ^ 

is a further coefficient to be evaluated numerically. This relation was missing in Paper I because the thermal equilibrium of 
the disc had not been considered. It may be helpful to point out that this relation, but with Q5 = 1, corresponds to the usual 
relation between the vertically integrated viscosity vYj and the surface density in a flat disc. The two are linked in that case 
by vYj — ctlQ.. One of the advantages of using T is that it allows a consideration of inviscid discs (see Paper I). The coefficient 
Qs reflects the thickening of the disc in response to increased viscous dissipation. 



3 STABILITY WITH RESPECT TO SHORT- WAVELENGTH PERTURBATIONS 

As in Paper I, the evolutionary equations derived are of the general form 

f + = «■ (69) 

+ ll^hl) = ~{^ + ^ d ^ +93 rl,^), (70) 

where v(r,t) is the mean radial velocity, h(r) = r 2 Q. is the specific angular momentum, and the torque components have the 
form gi(r, E, The system is of parabolic type and can be thought of as a non-linear advection-diffusion-dispersion system. 

In the case of a flat disc, the problem reduces to a non-linear diffusion equation for the surface density. It is well known 
(Lightman & Eardley 1974) that local thermal-viscous instability occurs if 9(z5£) /<9£ < 0. (It is assumed that b! > 0, because 
otherwise the orbital motion itself is unstable.) In the notation of the present paper, this criterion corresponds to dl/dT, < 0. 
It is important to find the generalization of this condition to the case of a warped disc. 

First eliminate v to obtain 
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at r or 



d£ dl 
gii + gir— + g 3 rl x — 
or or 



dgi 
dr 



.92 



h_ 

1? 



(71) 



Consider the stability of any given solution of this equation with respect to linear Eulerian perturbations (dE, 8£). Note that 

^i 5E+ H^i' (72) 

while 



2 dl d5l 



> | dr dr 

and £ • 51 = 0. The perturbed angular momentum equation has the form 



(73) 



h^-(SXl + Z5l) = 
ot r or 



5 9l £ + 9l 8£ + 5g 2 r 



dl 



dr 

fd5g, 5(g 2 M 2 ] 
\ dr r 



92r 



d5l 
dr 



di 



di 



+ 5g 3 rl x — + g 3 r5£ x — + g 3 rl x 



dr 

dgi _ g 2 \ip\ 2 
dr r 



dr 



dSl 
dr 



61 



(74) 



Now suppose the perturbations vary much more rapidly with radius and time than the unperturbed solution, so that they 
have the form 

exp I — i / ui dt + i / k dr 



(75) 

multiplied by more slowly varying functions. Here \uit\, \kr\ S> 1 and one may formally let k — » oo (although in reality the 
equations become invalid when k is comparable to 1/H). If one adopts the arbitrary scaling \5l\ — 0(k~ 1 ), then <5E = O(l), 
6\tp\ = O(l) and 5gi = O(l), while w = 0{k 2 ). 

Equation ( [74] ) should then be projected on to the orthonormal basis (£,m,n) defined by r(dl/dr) = \ip \ m and Ixm = n. 
Writing 51 — 5mm + 5nn, one obtains, at leading order in k, 



iuihr 5T, 



-iuihrTjSm = \k\ib\ 5g 2 — k 2 g 2 r 5m + k 2 g 3 r 5n — ik\ip\ —r- Sgi, 

rh 



-\LohrY,5n = 



-k 2 g 2 r 5n + ik\ip\ 5g 3 — k 2 g 3 r 5m, 



(76) 

(77) 
(78) 

while 5\ip\ — \kr 5m. The solvability condition yields a dispersion relation in the form of a cubic equation for Lojk 2 ; instability 
occurs if any of the roots has a positive imaginary part, since then short-wavelength disturbances grow exponentially in time. 
It has been shown that, under the assumptions of this paper, the torque coefficients have the form 

Qi = Qi(\<ip\)Ir 2 n 2 , 



with 

T = Q 5 (|V!)S p F(r) 



(79) 



(80) 



where the exponent p = 5/3 (Thomson) or 10/7 (Kramers), and F(r) is some function. Define the dimensionless growth rate 
s by 



— iuj = s—Qk 2 ; 

then the cubic dispersion relation corresponds to 

s - apQi -a(QiQ s )'/Q 5 

det ' 







= 0, 



(Q2 - aQ 1 )p\iP\ s + Q 2 + (Q2Q5)'\ip\/Q5-a(QiQ 5 y\ip\/Q5 -Qs 
pQs\1>\ Qa + iQaQsYM/Qs s + Q 2 J 

where a = h/(rh'), and a prime on a Q coefficient denotes differentiation with respect to 
If the disc is initially flat [\tp\ — 0), the dispersion relation has roots 



apQ\ 



-Q 2 ± iQ 3 



(81) 



(82) 



(83) 



For stability, one requires apQi < and Q 2 > 0. Note that p > (or, rather, ap > 0) is the usual criterion for thermal- viscous 
stability as discussed above. For a flat disc, the Q coefficients always satisfy Qi < and Q 2 > 0, except possibly for unphysical 
rotation laws. The condition Q 2 > means that the effective diffusion coefficient for a small warp introduced into the disc is 
positive. 

For an initially warped disc, the dispersion relation contains as a parameter and its roots must be determined 
numerically along with the Q coefficients. Instability occurs when a root exists with Re(s) > 0. 
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Table 1. Critical warp amplitude for local instability in a Kcplcrian disc. 
a \tp\ (Thomson) \ip\ (Kramers) 



0.01 0.19 0.17 

0.02 0.29 0.24 

0.03 0.37 0.31 

0.04 0.47 0.39 

0.05 0.63 0.47 



4 NUMERICAL RESULTS 

The numerical evaluation of the dimensionless torque coefficients proceeds similarly to Paper I. Attention is restricted here 
to the case of a Keplerian disc (g = 3/2) with adiabatic exponent T = 5/3 and vanishing bulk viscosity («b = 0). 

In Figs 1-3 the coefficients Qi, Q2 and Q a , and also the products Q1Q5, Q2Q5 and Q3Q5, are plotted as functions of 
the amplitude \tp\ of the warp. Three values for the viscosity coefficient, a = 0.01, 0.1 and 1, are investigated, and both 
Thomson and Kramers opacity are considered. The coefficients Qi, Q2 and Qz depend very little on the opacity law and are 
very similar to their values for the polytropic disc considered in Paper I. However, the coefficient Q 5 increases rapidly with 
increasing amplitude, reflecting the thickening of the disc in response to increased viscous dissipation. The precise behaviour 
of Q5 depends to some degree on the opacity law. The products Q1Q5, Q2Q5 and Q3Q5 reflect the true variation of the three 
torque components with the amplitude of the warp. 

If a is large, say a = 1, the variation of Q1Q5, Q2Q5 and Q3Q5 with amplitude is modest. However, for smaller a the 
variation can be significant. In particular, if a Ss 0.05, Qi can eventually become positive, meaning that the usual viscous 
torque parallel to I is reversed. However, as noted in Paper I, the disc becomes highly non-axisymmetric in this limit and the 
validity of the solution appears questionable. 

The roots of the local dispersion relation for short-wavelength perturbations can be calculated from this information. The 
system is found to be stable for < |V>| < 1 when a = 0.1 or a — 1. However, for smaller a, instability is found for sufficiently 
large \tp\. The critical values of \%p\ for various values of a are listed in Table 1. Instability usually, although not invariably, 
occurs where Q\ is already positive. 



5 SUMMARY AND DISCUSSION 
5.1 Summary of results 

Under the assumptions of this paper, the non-linear dynamics of a warped accretion disc may be described by one-dimensional 
conservation equations for mass, 

f + ^)=0, (84) 
and angular momentum, 

9 {Er 2 ne) + id 2 Q£ l^ +T ^ (g5) 
at r Or r Or 

Here £(r, t) is the surface density, v(r,t) the mean radial velocity, fi(r) the orbital angular velocity and £(r,t) the tilt vector. 
The internal torque is 2irG(r, t), given by 

G = Qilr 2 ti 2 I + Q 2 lr 3 S? ^ + Q 3 lr 3 n 2 fx^, (86) 

or or 

where T(r, t) is the azimuthally averaged second vertical moment of the density. Finally, T(r, i) represents any external torque 
due to self-gravitation, radiation forces, tidal forcing, etc. The dimensionless coefficients Qi, Q 2 and Q 3 depend on the rotation 
law and the shear viscosity, and also, in the non-linear theory, on the adiabatic exponent, the opacity law, the bulk viscosity 
and the amplitude |V>| = r\di/dr\ of the warp. 

For an ideal gas and a power-law opacity function, the relation between I and E is 

T= Q 8X . E (10+S-2»)/(6+«- 2w ) ( ^2 a)a /(6+«- a v) n -2(5-2 1( )/(6+ B -2v) ^m^H ^ -2(4-y)/(6+*-2„) -2/(6+x-2 y ) ^ ^ 

where Q5 is a further dimensionless coefficient, equal to unity for a flat disc. Rosseland mean opacities for a standard 'solar' 
chemical composition have been computed by, e.g., Iglesias & Rogers (1996), who tabulate log kh as a function of logT and 
log R = log p — 3 log T + 18 (all quantities expressed in CGS units). Power-law opacity functions of the Thomson and Kramers 
type offer a fair approximation to two large sectors of this table. Fiducial values of C K = 0.33 for the Thomson region (based 
on logT = 7 and logi? = —7) and C K = 4.5 x 10 24 for the Kramers region (based on logT = 4.5 and logi? = —3) will be 
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Figure 1. The variation of the coefficient Qi (dashed lines) and the product Q1Q5 (solid lines) with the amplitude \ip\ of the warp, for 
a Kcplcrian disc with T = 5/3 and no bulk viscosity. The left-hand panels are for Thomson opacity, the right-hand panels for Kramers 
opacity. The upper, middle and lower panels correspond to a = 0.01, 0.1 and 1 respectively. 
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Figure 2. The coefficient Q2 (dashed lines) and the product Q2Q5 (solid lines), as for Fig. 1. 
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Figure 3. The coefficient Q3 (dashed lines) and the product Q3Q5 (solid lines), as for Fig. 1. 
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adopted. The errors in these quantities, especially the latter, are quite large, but fortunately the dependence of 1 on C K is 
very weak (T oc C\J A for Thomson and X oc ClJ 7 for Kramers), so C K need not be defined very accurately. 
For a Keplerian disc (q = 3/2) with /u, m = 0.6, this gives 

I w 6.5 x 10 11 Qb a 1/3 E 5/3 fT 5/3 (88) 

for Thomson opacity and 

J ?s 4.4 x 10 12 Qs a 1/7 E 10/7 !^ 12/7 (89) 

for Kramers opacity, where all quantities are referred to CGS units. The numerical coefficients quoted here involve estimated 
uncertainties of 10-20% owing to variations in C K and /i m , depending on how the power-law opacity regions are delineated. 
Minor adjustments may also be required to account for the effects of finite optical thickness and/or external irradiation, as 
discussed in the Appendix. 

The Q coefficients can readily be calculated numerically by solving a set of ODEs in azimuth. These are 

f = (3-y)/5/i//2-(6 + a;-2y)/i/4, (90) 
fi = (r + l)/ 2 /4 + (r-l)/ 2 {2a(/3|^|cos^) 2 + 2a(/4 + /3|V|cos^.) 2 

+ a [(fi + h\i>\ cos0) \ip\ cos(j)- \ip\ sin0 - / 3 ] 2 + a (-q + f r ,\tp\ cos<^>) 2 + a/5 + (a b - §a) ft - g 2 a/i} , (91) 

f = hh + 2fr+ [l + (a b + I a )/ 4 ]/ 2 ^|cos (? ! ) -a/ 2 / 3 (l + |V| 2 cos 2 0)-Q/ 2 |V|sin0, (92) 
/I = -/3|^|cos0 + 2/3|i/»|sin0 + / 4 (/4 + /3|V|cos0) + l- [l + (a b + |a)/ 4 ] / 2 

—a/2 (/4 + fz\ip\ cos0) (1 + |i/>| 2 cos 2 (f>) + a/ 2 |i/>| 2 cos sin 0, (93) 

/s = /4/b - (2 -g)/a- 0/2/5(1 + H 2 cos 2 0) +ga/ 2 |V|cos^, (94) 

/e = -2/4/6, (95) 

in which /„ stands for fn{4>). These are subject to periodic boundary conditions /n(27r) = f n (0) and the normalization 
condition (fe) = 1, where the angle brackets denote azimuthal averaging. The required coefficients follow from 

Qi = (/6[-?a/2-/ 3 /5 + a/ 2 / s |V|cos0]), (96) 
<3 2 +iQ3 = |i|(e"*/ 6 [/3 - i/3(/4 + f:i\^\ cos 0) + ia/ 2 (/4 + /3M cos </>)|^| cos (/> - ia/ 2 /3 - ia/2 1^| sin0] ), (97) 

Q 5 = ( / 2/(6+*-2y) / -2(3- ! ,)/(6 +a; -2 ! ,)^ (gg) 

Apart from the torque component with coefficient Q3, a direct correspondence with the notation of Pringle (1992) can 
be made by identifying 

Q 1 1Q «-» (99) 
dlnr 

Q 2 Jf7 <-> ^ 2 E. (100) 

Therefore 1/1 oc —Q1Q5 and i^ 2 oc Q 2 Q5. The numerical evaluation of the Q coefficients (Section 4) demonstrates that the 
qualitative behaviour of v\ and i^ 2 depends on the value of the small-scale viscosity coefficient a. For large a, say a = 1, both 
v\ and ^ 2 increase with increasing amplitude. For smaller a, v\ decreases with increasing amplitude and may even become 
negative; while V2 first increases and then decreases. Also important is the variation of v\ and V2 with S, which is the same 
as for a flat disc. 



5.2 Discussion 

The non-linear theory derived in this paper is suitable for many applications. Some comment is warranted on the meaning of 
the solution and the conditions under which it is expected to be valid. 

Starting from the basic fluid-dynamical equations in three dimensions, I have derived one-dimensional evolutionary 
equations using a consistent asymptotic expansion that is based on the thinness of the disc, H/r = O(e), and the consequent 
separation of the orbital and evolutionary time-scales. The theory is fully non-linear in the sense that warp amplitudes \tjj\ of 
order unity are allowed for. Although this 'solution' is not an exact solution of the Navier-Stokes and accompanying equations, 
the theory presented here is asymptotically exact, meaning that it differs from an exact solution by an amount that tends to 
zero as O(e) in the limit e — ► 0. It may be noted that any attempt to define the tilt vector of a disc of non-zero thickness to 
an accuracy better than O(e) is somewhat arbitrary. 

The validity of the solution should be questioned when this error, although formally O(e), is numerically large and 
accumulates rapidly. This is expected to occur for a Keplerian disc when a <S H/r, when a resonance occurs, as described in 
Paper I. In this case e is not small enough for the asymptotic analysis to apply. On the other hand, the SPH simulations of 
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Nelson & Papaloizou (1999) suggest that, even when a ^ H/r, the evolution of warps again becomes diffusive in character 
when the amplitude is large. This is consistent with the observation made in Paper 1 that the resonant behaviour is weakened 
at large amplitude; it can also be seen from the upper panels of Fig. 2, where the resonant peak in the coefficient Q% is located 
at small amplitude only. It would be useful to make more detailed comparisons between SPH simulations and the present 
theory. 

Although Nelson & Papaloizou (1999) found shock-like features in some of their simulations, this occurred when the 
wavelength of the bending disturbance was comparable to the disc thickness so that a collision occurred between oppositely 
directed horizontal flows. The present theory does not allow for such short-wavelength disturbances. In any case, it is not 
expected that strict discontinuities will occur in a solution of the Navier-Stokes equation. 

A further issue is the stability of the solution. As shown in Section 3, if a fiat disc is stable to the thermal-viscous 
instability, it usually remains viscously stable when it is warped.^ An exception arises if a is small, say a <J 0.05, when the 
disc may become locally unstable at sufficiently large warp amplitude. It has been noted, however, that the deformation of 
the disc is more extreme in this limit, and the induced motions are more likely to be hydrodynamically unstable. Therefore 
the viscous instability, whose physical origin is unclear, may not occur in practice. It is not evident what the true outcome 
of a large-amplitude warp in a low-viscosity disc would be, but the parametric instability discussed by Gammie, Goodman & 
Ogilvie (2000) is likely to be important. 

Provided that a is not so small as this, the principal remaining uncertainty in this analysis concerns the modelling of 
the turbulent stress tensor. The simplest assumption, adopted here as in almost all theoretical work on accretion discs (and, 
implicitly, in SPH simulations), is that the stress behaves similarly to an isotropic viscosity. Measurements of the damping 
rate, caused by magnetohydrodynamic turbulence, of shearing motions of the type induced in a warped disc (Torkelsson et al. 
2000) broadly support this hypothesis, but more numerical and theoretical effort is required in this area. 

The very innermost part of a disc around a neutron star or black hole is dominated by radiation pressure and Thomson 
scattering opacity. This case would require a treatment of the radiation-hydrodynamic equations and is beyond the scope of 
this paper. However, one can anticipate that the evolutionary equations are identical in form and that the appropriate X-E 
relation for a flat disc is modified with a Q5 coefficient in a similar way. However, the detailed variation of the Q coefficients 
would require a further investigation. 

Discs in binary systems are subject to tidal forces that cause a distortion of the flow, possibly including spiral shocks. The 
interaction of such distortions with a warp requires further investigation. If spiral shocks are the dominant mode of angular 
momentum transport in the disc, the 'viscous' theory developed here and elsewhere may not be relevant. 

In conclusion, this paper provides the logical generalization of the alpha theory of Shakura & Sunyaev (1973) to the case 
of a time-dependent warped accretion disc. It is derived directly from the basic fluid-dynamical equations, and is valid even 
for warps of large amplitude. The scheme is suitable for practical applications: it is expressed in physical units, it requires the 
solution of only one-dimensional equations from which the fast orbital time-scale has been eliminated, and it allows numerous 
external effects to be included in terms of the torque that they exert on the disc. Such applications will be presented in future 
publications. 
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APPENDIX A: EFFECTS OF FINITE OPTICAL THICKNESS AND EXTERNAL IRRADIATION 

In Sections 2.3 and 2.4 the vertical structure of the disc was determined by applying the 'zero boundary conditions' p — T — 
at the surface. This is correct in the limit of a highly optically thick disc. In practice the optical thickness is usually large 
but certainly finite, and this leads to a small correction to the dynamical equations through a change in the dimensionless 
quantity X* . More importantly, the disc may be subject to an external irradiating flux, which can affect the vertical structure. 
This is particularly important for discs around neutron stars and black holes, ff the disc is warped it is likely to intercept a 
significant fraction of the radiation emitted by the central X-ray source. 

The effects of finite optical thickness and external irradiation may be estimated from a simple model, derived in part from 
a more detailed analysis by Dubus et al. (1999). It is supposed that the equations of Section 2.2 extend up to a photospheric 
surface C = Csj above which F = F s — constant and T — T s = constant. The constancy of F implies that there is negligible 
dissipation of energy in the atmosphere, while the constancy of T is a good approximation to more detailed atmospheric 
models. The possibility of a hot corona is neglected here. 

The photospheric boundary conditions are taken to be 

F s + = aT B \ (Al) 

and 

rs = §, (A2) 
where Ti rr is the irradiation temperature and 

KRprdC (A3) 

the optical depth of the photosphere. To a good approximation this equates to 

'C-C 







d(C-Cs), (A4) 



where h s is the density scale-height (in terms of the variable £) at the photosphere, given by 
kT s 1 / 1 dp' 



/t* m TOH h s \ pd( / 
This may be rearranged into the form 



(A5) 



For a fiat disc, the analysis of Section 2.3 proceeds as before, but now with boundary conditions 

<SF,(C s ,)+ITL = [T»(Cs,)] 4 (A7) 
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Table Al. Effects of finite optical thickness and external irradiation on the second moment in the case of a fiat disc with Kramers 
opacity. 



X 



71 
irr* 


T 1 

-* s* 


7 1 

-t c* 




T 


0.001 





0.137 


0.802 


2740 


0.708 


0.01 





0.244 


0.803 


274 


0.699 


0.1 





0.437 


0.814 


27.0 


0.622 


0.001 


0.3 


0.303 


0.803 


2710 


0.704 


0.01 


0.3 


0.329 


0.804 


271 


0.691 


0.1 


0.3 


0.460 


0.816 


26.8 


0.614 


0.001 


1.0 


1.000 


1.027 


723 


0.899 


0.01 


1.0 


1.001 


1.031 


77.7 


0.753 


0.1 


1.0 


1.013 


1.049 


9.61 


0.503 



and 

[p,(C B *)] 1+x p;(Cs*)] 1+H = ±«S(1 + *)£., (A8) 
where 

T irr = r irr „ Ut, (A9) 



and 

c Uf 



all* 



(A10) 



is a small dimensionless parameter. When the solution is obtained, the optical depth at the mid-plane can be determined from 



2 f Cs 



l + Xrpy 



r, ^ / ,, riP rdC=| + g / /> M "T : "<lc (All) 



o 



Therefore 5 is an approximate inverse measure of r c . 

The computed values of X* for some representative values of S and Ti rr » are given in Table Al. Also given are the optical 
depth at the mid-plane r c , the surface temperature T s * and the central temperature T c *. The contribution of the atmosphere 
to £ and X has been neglected. It can be seen that a noticeable change in X* occurs, as expected, when the optical thickness 
is not large, or when the irradiation temperature exceeds the central temperature of the non-irradiated solution. However, 
the fractional change is less than 30% even for a strongly irradiated disc with r c ~ 10. 

It has not been found possible to give an exact solution for the vertical structure of a warped disc including the photospheric 
boundary conditions. The difficulties are especially pronounced in the case of a disc irradiated by the central X-ray source, 
since the irradiation then varies strongly with azimuth (especially if the disc is partially self-shadowed) and is one-sided, 
breaking the symmetry of the vertical structure. However, for a < 1 the thermal time-scale is longer than the orbital time- 
scale and it is appropriate to average the irradiation over the orbital motion of the fluid. This restores the symmetry about 
the mid-plane and suggests that the effect on the T-E relation will be similar to that found above for a flat disc. 
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